%% STATIC: Optimally design MPAs wrt profit-biomass tradeoff - % MPA NOT
%% SPECIFIED
disp('STATIC1 = Everything static.')
%Values:
%Biomass: virgin biomass (Bj_initial) in MPAs + biomass in fished patches with optimal Fsum and no MPAs (Bi_Fsum_opt) in fished patches
%Profit: profit with optimal Fsum and no MPAs (Profit_i_Fsumopt_NoMPAs) in fished patches

%scaler for tracing EF (0=all conservation; 1=all payoff)
% Payoff_Conservation_EF_alpha_STATIC1_range=[0:0.2:1]; %test run
% Payoff_Conservation_EF_alpha_STATIC1_range=[0:0.1:1]; %lo res
% Payoff_Conservation_EF_alpha_STATIC1_range=linspace(0,1,15);
Payoff_Conservation_EF_alpha_STATIC1_range=linspace(0,1,25); %added more values (was 15 above)
Payoff_Conservation_EF_alpha_STATIC1_range(end)=0.99; %fix1
% Payoff_Conservation_EF_alpha_STATIC1_range=[0 0.125 0.25:0.05:0.55 0.7 0.8 0.9 0.99];% lower N
% Payoff_Conservation_EF_alpha_STATIC1_range=[0:0.05:1]; %med res
% Payoff_Conservation_EF_alpha_STATIC1_range=[0:0.01:1]; %publication

%shells
Payoff_Conservation_EF_wrt_alpha_STATIC1_naive=NaN(2,length(Payoff_Conservation_EF_alpha_STATIC1_range)); %storing payoff and conservation outcomes
Payoff_Conservation_EF_wrt_alpha_STATIC1_informed=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive;
tictoc_wrt_Payoff_Conservation_EF_wrt_alpha_STATIC1=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range)); %time keeper
MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1=NaN(length(domain_suit_hab),length(Payoff_Conservation_EF_alpha_STATIC1_range)); %for storing MPA locs solns
BestObjective_wrt_alpha_STATIC1=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range));
PointHistory_wrt_alpha_STATIC1=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range));
PointObjectiveHistory_wrt_alpha_STATIC1=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range));
func_evals_wrt_alpha_STATIC1=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range));
Fi_EF_wrt_alpha_STATIC1=NaN(length(domain_suit_hab),length(Payoff_Conservation_EF_alpha_STATIC1_range));

%fill shell with extreme (alpha=0 [conservation] and 1 [profit])
if NPVon==0 %equil
    %note: assuming that Fsum opt produces optimal MPA design for maximizing profit (may be 0% MPAs)
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,1)=0; %payoff with alpha=0
    % Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,1)=Kvirgin; %conservation with alpha=0 %NO v7
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,1)=sum(Kunfished_hat); %conservation with alpha=0 %what naive manager thinks will happen based on biomass in existing MPAs v7
    % Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,end)=FleetProfit_Fsum_opt; %payoff with alpha=1
    % Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,end)=Bijy_Fsum_opt; %conservation with alpha=1
    % Payoff_Conservation_EF_wrt_alpha_STATIC1_informed=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive; %for alpha=0 naive=informed %NO v7
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(1,1)=0; %payoff will be zero (naive manager got this correct) v7
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(2,1)=Kvirgin; %actual biomass that will happen if closed entire domain into MPAs v7
    MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1(:,1)=1; %MPA design with alpha=0
    % MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1(:,end)=(Fi_Fsum_opt==0); %MPA design with alpha=1
    Fi_EF_wrt_alpha_STATIC1(:,1)=0; %no harvest bc all MPAs
    % v7:
    %for alpha=1 (all profit), leave system as is (no new MPAs, because the
    %static manager invariably considers MPAs a cost)
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,end)=sum(Profit_i_Fsumopt_NoMPAs); %existing profit with no new MPAs
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,end)=sum(Bi_Fsum_opt); %existing biomass at Fsum msy
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(:,end)=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(:,end); %informed result the same, bc nothing changed
    MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1(:,end)=NCI_MPA_locs';%existing MPAs CORRECTION 1/7/14
    EXITFLAG_static=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range));
elseif NPVon==1
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,1)=0; %payoff with alpha=0
    % discount Kunfished_hat
%     Kunfished_hat_discount_iy=repmat(Kunfished_hat,1,max_simulation_length).*discounted_value_iy;
%     Kunfished_hat_discount=sum(Kunfished_hat_discount_iy(:));
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,1)=Kunfished_hat_discount; %conservation with alpha=0 %what naive manager thinks will happen based on biomass in existing MPAs v7
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,1)=sum(Kunfished_hat);
    %informed with all MPAs (ie Fsum=0)
    Fsum=0;
    MOPmodel2
    tmp1=sum(Profit_iy,1);
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(1,1) =tmp1(end); %payoff
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(2,1)=sum(sum(Bijy(:,:,end))); %conservation; %actual biomass that will happen if closed entire domain into MPAs v7
    MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1(:,1)=1; %MPA design with alpha=0
    Fi_EF_wrt_alpha_STATIC1(:,1)=0; %no harvest bc all MPAs
    %commented out below (fix1)
%     %for alpha=1 (all profit), leave system as is (no new MPAs, because the
%     %static manager invariably considers MPAs a cost)
%     Fsum=sum(Fi_Fsum_opt);
%     MPA_locs=NCI_MPA_locs;
%     MOPmodel2
%     tmp1=sum(Profit_iy,1);
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,end)=tmp1(end); %payoff
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,end)=sum(sum(Bijy(:,:,end)));
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(:,end)=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(:,end); %informed result the same, bc nothing changed
%     MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1(:,end)=NCI_MPA_locs';%existing MPAs CORRECTION 1/7/14
%     EXITFLAG_static=NaN(1,length(Payoff_Conservation_EF_alpha_STATIC1_range));
    
end

%now find payoff-conservation EF for every 0<alpha<1  
% x0_Fsumguess=Fsum_opt; %use as starter, but then update wrt alpha and MPAlocs
Payoff_Conservation_EF_alpha_STATIC1_i=1;
tic

%heuristic:
Dims=numpatches-length(existMPAs); % #patches not including those fixed as MPAs
MaxSingleChangeTrials=1e4; %default 1e4
TabuOnOff=1; %default=1
%Crow: have it start at reasonable #MPAs:
Point = zeros(Dims,1); % Initiate point
idx = randperm(Dims); idx = idx( 1:(floor((1-Payoff_Conservation_EF_alpha_STATIC1_range(2))*Dims)));
Point( idx ) = 1; Point_initial=Point;

% for Payoff_Conservation_EF_alpha_STATIC1=Payoff_Conservation_EF_alpha_STATIC1_range(2:end);
% for Payoff_Conservation_EF_alpha_STATIC1=Payoff_Conservation_EF_alpha_STATIC1_range(2:(end-1)); %v7
for Payoff_Conservation_EF_alpha_STATIC1=Payoff_Conservation_EF_alpha_STATIC1_range(2:end); %fix1
    Payoff_Conservation_EF_alpha_STATIC1_i=Payoff_Conservation_EF_alpha_STATIC1_i+1;
    
    disp(['--Payoff_Conservation_EF_alpha_STATIC1 = ',num2str(Payoff_Conservation_EF_alpha_STATIC1),])
%     Fsum=Fsum_opt; %that is if do not optimize Fsum each time
       
    t1start = tic;
%     [ BestPoint, BestObjective, PointHistory, PointObjectiveHistory, func_evals ]=...
%         BinaryTabuMonteCarlo(@Payoff_Conservation_EFvalue_STATIC1, Dims, MaxSingleChangeTrials,TabuOnOff, Point_initial);

% genetic algorithm
%  [X,FVAL,EXITFLAG,OUTPUT]=ga(@Payoff_Conservation_EFvalue_STATIC1,Dims,[],[],[],[],0,1,[],ones(1,Dims));
% options_ga=optimset('Display','off','InitialPopulation', Point_initial');
options_ga=gaoptimset('Display','off','InitialPopulation', Point_initial');
 [X,FVAL,EXITFLAG,OUTPUT]=ga(@Payoff_Conservation_EFvalue_STATIC1,Dims,[],[],[],[],zeros(1,Dims),ones(1,Dims),[],[1:Dims],options_ga);
 BestPoint=X';
 BestObjective=FVAL;
 func_evals=OUTPUT.funccount;
 

%%
    Point_initial=BestPoint; %update heuristic starting point
    BestObjective_wrt_alpha_STATIC1(Payoff_Conservation_EF_alpha_STATIC1_i)=BestObjective;
%     PointHistory_wrt_alpha_STATIC1(Payoff_Conservation_EF_alpha_STATIC1_i)=PointHistory;
%     PointObjectiveHistory_wrt_alpha_STATIC1(Payoff_Conservation_EF_alpha_STATIC1_i)=PointObjectiveHistory;
    func_evals_wrt_alpha_STATIC1(Payoff_Conservation_EF_alpha_STATIC1_i)=func_evals;
    EXITFLAG_static(Payoff_Conservation_EF_alpha_STATIC1_i)=EXITFLAG;

    tmp=NCI_MPA_locs'; %set vector just with NCI MPAs
    tmp(NCI_MPA_locs==0)=BestPoint; %fill in the new MPAs by heuristic
    BestPoint_allMPA_locs=tmp;
    MPAlocs_Payoff_Conservation_EF_wrt_alpha_STATIC1(:,Payoff_Conservation_EF_alpha_STATIC1_i)=BestPoint_allMPA_locs; %store MPA locs soln
      
    %naive static results (expected, though not really going to happen)
    MPA_locs_STATIC1=BestPoint_allMPA_locs; MPA_locs=BestPoint_allMPA_locs'; %locs that are fished (0) or MPAs (1) Q: NEED TO TRANSFORM
%     Bi_mpas=sum(Bi_initial.*MPA_locs_STATIC1); %biomass in the MPAS (=virgin biomass in those patches) - Manager does not know this
    Bi_mpas=sum(Kunfished_hat.*MPA_locs_STATIC1); %biomass in the MPAS (=regtression of NCI MPAs interpolation in those patches)
   
    Bi_fished=sum(Bi_Fsum_opt.*(MPA_locs_STATIC1==0)); %biomass in fished patches (=biomass in those patches given optimal Fsum and no MPAs)
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,Payoff_Conservation_EF_alpha_STATIC1_i)=Bi_mpas+Bi_fished;
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,Payoff_Conservation_EF_alpha_STATIC1_i)=sum(Profit_i_Fsumopt_NoMPAs.*(MPA_locs_STATIC1==0)); %payoff (=profit in fished patches=profit in those patches given optimal Fsum and no MPAs)
    
    %informed static results
    %manager will set Fsum optimally in response to chosen MPAs and
    %dynamics (biophysical and fleet model)
    if length(find((MPA_locs==1)))==numpatches %    100% MPAs
        Fsum=0;
        Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(1,Payoff_Conservation_EF_alpha_STATIC1_i)=0; %payoff
        Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(2,Payoff_Conservation_EF_alpha_STATIC1_i)=Kvirgin; %conservation
    else        
    %npsol==================================================
        TolX_value=1e-4*numpatches; %Set tolerance for Fsum search. 1e-4 is the default TolX. 
    TolFun_value=TolX_value;
    options_fmincon_OptFsum=optimset('Algorithm','interior-point','Display','off','TolX',TolX_value,'TolFun',TolFun_value);
%     options_fmincon_OptFsum=optimset('Algorithm','interior-point','Display','off');
    prob.obj = 'OptFsum';
    prob.lo = [0];
    prob.up = [Fsum_UB_prop_wrt_sumFishingHab*sum(domain_suit_hab(MPA_locs==0))];
    prob.guess=(prob.lo+prob.up)/2;
%     [NEG_FleetProfit_Fsum_opt,Fsum_opt, g, hess, flag, iter]=npsol(prob);  % OPTIMIZATION
    [Fsum_opt,NEG_FleetProfit_Fsum_opt,flag,output]=fmincon(@OptFsum,prob.guess,[],[],[],[],prob.lo,prob.up,[],options_fmincon_OptFsum); %v7
    if Fsum_opt>= prob.up && sum(domain_suit_hab(MPA_locs==0))>0; disp('WARNING: Fsum_opt>= prob.up > 0'); end
    %=========================================================
    Fsum=Fsum_opt;
    MOPmodel2
    tmp1=sum(Profit_iy,1);
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(1,Payoff_Conservation_EF_alpha_STATIC1_i)=tmp1(end); %payoff
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(2,Payoff_Conservation_EF_alpha_STATIC1_i)=sum(sum(Bijy(:,:,end))); %conservation
    Fi_EF_wrt_alpha_STATIC1(:,Payoff_Conservation_EF_alpha_STATIC1_i)=Fi;
%     data_save_midrun=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed;
%     save data_save_midrun
    end
    tictoc_wrt_Payoff_Conservation_EF_wrt_alpha_STATIC1(Payoff_Conservation_EF_alpha_STATIC1_i)=toc(t1start)/60;
    disp(['#MPAs=',num2str(length(find(MPA_locs==1))),'; Fsum=',num2str(Fsum),'; took ',num2str(tictoc_wrt_Payoff_Conservation_EF_wrt_alpha_STATIC1(Payoff_Conservation_EF_alpha_STATIC1_i)),' minutes; EXITFLAG=',num2str(EXITFLAG)])
    disp(['Payoff_Conservation_EF_alpha_STATIC1=',num2str(Payoff_Conservation_EF_alpha_STATIC1),'; out (goal is NEG ONE)=',num2str(BestObjective), ]);
end
disp(['Minutes of each heuristic: ',num2str(tictoc_wrt_Payoff_Conservation_EF_wrt_alpha_STATIC1),])
disp(['Total hours of entire EF heuristic: ',num2str(nansum(tictoc_wrt_Payoff_Conservation_EF_wrt_alpha_STATIC1)/60),])
%%
figure %absolute outputs
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(2,:),Payoff_Conservation_EF_wrt_alpha_STATIC1_naive(1,:),'ko','markerfacecolor','k')
hold on
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(2,:),Payoff_Conservation_EF_wrt_alpha_STATIC1_informed(1,:),'r^','markerfacecolor','r')
xlabel('Biomass [kg]')
ylabel('Profit [$]')
set(gcf,'color',[1 1 1]); box off
legend('Naive static','Informed static')

%scale outputs
% if NPVon==0 %equil
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed;
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:)./FleetProfit_Fsum_opt; %scale payoff to fleet model max (note: not necessarily 0% MPA)
    Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:)./Kvirgin; %scale conseration to 100% MPA outputs
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive;
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:)./FleetProfit_Fsum_opt; %scale payoff to fleet model max (note: not necessarily 0% MPA)
    Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:)./Kvirgin; %scale conseration to 100% MPA outputs
% elseif NPVon==1 %npv
%     
%     Fsum=sum(Fi_Fsum_opt);
%     MPA_locs=NCI_MPA_locs;
%     MOPmodel2
%     tmp1=sum(Profit_iy,1);
%     FleetProfit_Fsum_opt=tmp1(end); %payoff
%     
%     Fsum=0;
%     MOPmodel2
%     Kvirgin=sum(sum(Bijy(:,:,end)));
%     
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed;
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:)./FleetProfit_Fsum_opt; %scale payoff to fleet model max (note: not necessarily 0% MPA)
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:)./Kvirgin; %scale conseration to 100% MPA outputs
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive;
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:)./FleetProfit_Fsum_opt; %scale payoff to fleet model max (note: not necessarily 0% MPA)
%     Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:)=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:)./Kvirgin; %scale conseration to 100% MPA outputs
% end

figure %scaled outputs
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:).*100,'ko','markerfacecolor','k')
hold on
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:).*100,'r^','markerfacecolor','r')
xlabel('Biomass [%]')
ylabel('Profit [%]')
set(gcf,'color',[1 1 1]); box off
legend('Naive static','Informed static')
ylim([0 110])
% %indicate alpha next to point
% for tmp=1:length(Payoff_Conservation_EF_alpha_STATIC1_range)
%     text_offset_x=.5;
%     text_offset_y=1;
%     text(Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,tmp).*100+text_offset_x,Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,tmp).*100+text_offset_y,num2str(Payoff_Conservation_EF_alpha_STATIC1_range(tmp)))
% end

%lines connecting naive and informed results
figure %scaled outputs
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:).*100,'ko','markerfacecolor','k')
hold on
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:).*100,'r^','markerfacecolor','r')
plot([Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:).*100;Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:).*100],...
    [Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:).*100;Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:).*100],'b')
hold on
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled(1,:).*100,'ko','markerfacecolor','k')
plot(Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(2,:).*100,Payoff_Conservation_EF_wrt_alpha_STATIC1_informed_scaled(1,:).*100,'r^','markerfacecolor','r')
xlabel('Biomass [%]')
ylabel('Profit [%]')
set(gcf,'color',[1 1 1]); box off
legend('Naive static','Informed static')
ylim([0 110])


data_save=Payoff_Conservation_EF_wrt_alpha_STATIC1_naive_scaled;
save data_save
